Supersolid phase transitions for hardcore bosons on a triangular lattice 
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Hard-core bosons on a triangular lattice with nearest neighbor repulsion are a prototypical ex- 
ample of a system with supersolid behavior on a lattice. We show that in this model the physical 
origin of the supersolid phase can be understood quantitatively and analytically by constructing 
quasiparticle excitations of defects that are moving on an ordered background. The location of the 
solid to supersolid phase transition line is predicted from the effective model for both positive and 
negative (frustrated) hopping parameters. For positive hopping parameters the calculations agree 
very accurately with numerical quantum Monte Carlo simulations. The numerical results indicate 
that the supersolid to superfluid transition is first order. 
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I. INTRODUCTION 

A supersolid is characterized by a non-trivial solid or- 
der parameter which co-exists vifith an ofF-diagonal long- 
range superfluid order parameter;^ This state appears 
counter-intuitive since the superfluid density implies bal- 
listic coherent transport without disturbing the ordered 
background of the solid. The possibility of a supersolid 
was discussed in the early 1970's and since then much 
attention has been focused on *He as the most likely ex- 
perimental realization, which remains controversial.^ The 
motions of many different kinds of defects have been con- 
sidered as possible mechanisms for supersolidity in ^He, 
including vacancies,^ grain boundaries'^ or dislocations^ 
which may be able to move through a solid state. In 
recent years supersolid behavior on lattices has been es- 
tablished theoretically in a number of interacting boson 
models, ^^^^ which now gives hope for an experimental 
realization of a supersolid with the help of ultra-cold 
gasesji2 One prototypical example to exhibit a supersolid 
phase are hard-core bosons with nearest neighbor repul- 
sion on a triangular lattice^"— which are described by 
the Hamiltonian 



^0 = -tj2{bib, 



/i^ni,(l) 



where foj (bi) is the boson creation (annihilation) opera- 
tor, t describes the hopping amplitude, /i is the chemical 
potential, V reflects the repulsion effect, and (i, j) repre- 
sents nearest-neighbor sites. 

For positive hopping i > a supersolid phase has 
been shown using numerical simulations^^— and varia- 
tional wavefunctionsjSiS For negative hopping t < it 
is not possible to use quantum Monte Carlo (QMC) sim- 
ulations, but at least for the special case of half filling 
a supersolid phase has been confirmedi^°di While the 
existence of the supersolid phase has been firmly estab- 
lished in this model, it would be desirable to obtain a 
better understanding of the physical mechanism leading 
to the supersolid phase and an analytical description of 
the quantum phase transition. 
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FIG. 1: The phase diagram of hard-core bosons on the trian- 
gular lattice. The dashed lines (blue) are numerical QMC re- 
sults while the solid lines (red) are the analytical results from 
Eqs. ((3]) and Q. The arrows indicate the order in terms of the 
classical directions in the spin-wave theory. Inset: Overview 
of the entire particle-hole symmetric region. 



The phase diagram of the model is shown in Fig. [T] 
For < ^ < 3V and t = the solid order is described by 
a 1/3-ordered state, where the unoccupied sites form a 
honeycomb sublattice around the filled sites as indicated 
in Fig. [21 Due to particle-hole symmetry the equiva- 
lent behavior can be seen for fi/V > 3 in the form of 
a 2/3 filled state (see inset of Fig. [T]). The goal of this 
paper is an analytic quantitative description of the tran- 
sition to the supersolid phase which occurs for i 7^ 
around fj,/V = 3. The solid order is three-fold degen- 
erate according to the choice of sublattice, however this 
degeneracy is not important for the formation of the su- 
persolid state so we will only consider one representative 
1/3-filled state in what follows. Finite hopping t ^ 
causes quantum fluctuations in form of virtual excita- 
tions from the occupied sites onto the unoccupied ones 
and back. Therefore the ordered state is renormalized, 
e.g. the density of the unoccupied honeycomb lattice is 



now finite ^ 3P/AV^ to lowest order and the occupied 
density is reduced ^ 1 ~ 3t^ /2V^ . However, the delo- 
cahzation of these virtual excitations is not the driving 
mechanism for the transition into the supersolid state. 
Instead we conjecture that a different type of excitation 
starts to appear at the phase transition, namely " defects" 
consisting of additional occupied sites. These defects can 
only move on the honeycomb lattice of unoccupied sites, 
while the solid order remains surprisingly stable, which 
leads to an effective model. 



II. EFFECTIVE MODEL 

The corresponding effective Hamiltonian for hard-core 
boson defects on the honeycomb lattice can be derived 
in the strong coupling expansion. Ordinary perturbation 
theory has a diverging number of contributions, but it 
is possible to analyze the movement of a single defect 
by only taking those terms into account which actually 
involve the defect This results in the following effective 
Hamiltonian on the honeycomb sublattice 
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FIG. 2: Schematic description of the effective model on 
the unoccupied honeycomb sublattice (lines) where the defect 
(large central sphere) can hop to neighboring sites. Smaller 
spheres (blue) are occupied sites, which only participate in 
virtual excitations. 



solid and supersolid 



i.d 



(ij) 



(2) 

where rrii = a|aj and the effective hopping parame- 
ters are given by ii — t + t'^ /V for nearest neighbors, 
t2 = t'^/V for next nearest neighbors, and ^3 = t'^ /V for 
next-next nearest neighbors on the honeycomb sublattice 
as shown in Fig. [21 The second-order terms arise from 
a two-step hopping process, where a boson from the oc- 
cupied sublattice jumps on an unoccupied site and then 
the defect takes the place of the original boson, thereby 
effectively transporting the defect. Likewise, the effective 
chemical potential is changed /i = /i — 3V + 5P/2V since 
virtual excitations from the unoccupied sublattice and 
back are now affected by the defect, which results in an 
effective self-energy. Since the defects can move ballisti- 
cally and coherently on the unoccupied honeycomb sub- 
lattice according to the effective model, a finite density 
of such defects also explains the coexistence of a finite su- 
perfluid density on top of a relatively stable renormalized 
ordered state, which are in fact the central characteristics 
of the supersolid phase. 

The effective model in Eq. ([2]) is of course still an inter- 
acting hard-core boson model, but is much easier to treat 
since we only have to consider small densities of defects 
as will be seen below. In particular, according to our 
hypothesis the phase transition line is simply determined 
by the point when it is energetically favorable to create 
one single defect. For positive hopping t > the energy 
of a single defect is minimized in the uniform k = state 
with a total energy of = —3ti — 6^2 — 3^3 — fl, which 
must change sign at the phase transition line. By taking 
also third order processes into account this results in the 
following prediction for the phase transition line between 



W -3t 



2V 41/2 



(3) 



As shown in Fig. [T] this result agrees remarkably well 
with our QMC simulations, where we used the stochastic 
series expansion with parallel tempering.~ The errorbars 
for the data from the QMC simulations are smaller than 
the size of the symbols in all figures. For fi > 3V the 
particle hole symmetric line is given by /i+ — 6V — 11+ . 

A similar analysis can be made for negative hopping 
t < 0. Since the honeycomb lattice is bi-partite the 
positive nearest neighbor hopping can be restored by a 
simple Gauge transformation of —1 on one sublattice, 
but the next-nearest neighbor coupling remains negative. 
Therefore, the ground state energy of one defect is now 
E- = 3^1 — 6t2 + 3t3 — fl, and accordingly we obtain a 
phase transition line of the form 



3V + 3t 



5t^ 



2V 



(4) 



where we have again taken third order processes into ac- 
count. For negative hopping it is more difficult to treat 
the problem with numerical methods, so there are no 
independent checks for this prediction so far. However, 
considering the good agreement for the positive hopping 
case above, it is reasonable to expect that the same cal- 
culation also applies for negative hopping with similar 
accuracy. We therefore find that the supersolid region is 
stable but slightly narrower for negative hopping, which 
appears to be caused by the additional kinetic frustra- 
tion. 
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FIG. 3: The density on the sites with lower fiUing on the 
triangular lattice near the phase transition (symbols) for t — 
0.03V = 300), t = 0.05V (V/3 = 150) and t = Q.IV 
{VP = 100). The solid lines are the analytical prediction 
from Eq. 



III. ANALYTICAL AND NUMERICAL 
RESULTS NEAR THE PHASE TRANSITION 

In order to calculate the detailed behavior near the 
phase transition analytically it is now possible to use 
linear spin-wave theory on the effective Hamiltonian in 
Eq. ([2]). It should be noted that the system is at low 
filling (i.e. nearly saturated in terms of spins) and that 
the honeycomb lattice has two nonequivalent sites, but 
otherwise the spin- wave calculation follows the standard 
schemciii In particular, the hard-core bosons can be 
mapped exactly onto spins, which are expressed in a ro- 
tated frame according to the classical alignment. The 
classical angle relative to the z-direction is given by 
cose = B/{W/2 - J2did) in terms of the field B ^ 
fi — ZV . The spins are then re-expressed in terms of 
regular bosons, which has the advantage that the inter- 
action can be solved by a Bogoluibov transformation to 
linear order. This procedure is especially reliable in the 
low density limit near the phase transition that we are 
interested in. The details of this calculation can be found 
in the appendix. The density {p) = \ — {S^) on the hon- 
eycomb sublattice at T = is then given in terms of the 
Bogoluibov rotation angles 



{p) = ^-cos( 



2 n0 
k 



, (5) 



where tanh2^^^ 

k ■ 
2 , 



and 
with 

\lk\^ 



tanh 29 
A = sin^ e {2t 

3 + 4 cos 



lk\A/ it+\jk\{A~t)) 
7fe|A/ {t-\^k\{A-t)) 
+ V) /4 and 



3fcx 
2 



^/3ky 



2 COS (y/Sky) 



1/2 



In Fig. [3] we show the density on the unoccupied sub- 
lattice from QMC simulations on a triangular lattice with 
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FIG. 4: The superfluid density on the triangular lattice near 
the phase transition (symbols). The solid lines are the ana- 
lytical prediction from Eq. (|6}. 



24^ sites near the solid to supersolid phase transition at 
different hoppings (symbols). This can be compared to 
the spin-wave results for the density on the honeycomb 
lattice from Eq. ([5]) (lines), where the small background 
contribution 3t^/4V^ from virtual excitations was added. 
The virtual excitations are localized so this finite back- 
ground density does not contribute to the superfiuid den- 
sity which is also confirmed by the data in Fig. H) Here 
and in what follows the spin-wave results always include 
the additional third order shift in p, by llfi/AV'^ from 
Eq. (|3l). The good agreement not only confirms the va- 
lidity of the effective model on the honeycomb lattice in 
Eq. (I2|), but also provides an analytical tool to calcu- 
late the general behavior and correlations near the phase 
transition for both positive and negative hopping t. As a 
check we have also computed densities and correlations 
on the honeycomb lattice, which agree equally well with 
the original triangular model (not shown). The corre- 
lations between the occupied sublattice and the unoc- 
cupied honeycomb sublattice can be calculated by the 
modified perturbation theory discussed above, which re- 
sults in positive correlations between sublattices 
for i > 0, while the correlations cancel to second order 
for t < 0, i.e. the occupied sublattice always gives only a 
very weak contribution to the superfluid density as indi- 
cated by the spin directions in Fig.m 

The superfluid density on the unoccupied honeycomb 
sublattice can be calculated in the effective spin model 
as the response to a small twist i> in the x-y plane 
Ps = (P{H{v)) /(Pv. Expanding the Hamiltonian leads 
to H{y) = H + Y.,A {"dJI + t^d/^Td), where i^d is the 
projection of the twist along the corresponding bond di- 
rection, JJ = {i/2)td {StS'^d - ^t+d^r) is the spin cur- 
rent, and Td = {l/2)td (StS'^d + ^t^d^r) is the spin- 
kinetic energyJ^iiS In linear spin-wave theory the re- 
sponse to the spin current does not contribute in the low 
density limit and the superfluid is therefore approximated 
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by the correlation functions Td of neighbors connected by 
id 



Ps - ^{Ti + iT2 + 3T3) 



(6) 



which can be calculated in the rotated spin frame (see 
appendix). In Fig.|4]we compare the results from Eq. (|6]) 
on the honeycomb lattice with the superfluid density ps 
derived from the winding number {W"^) /Apt in QMC sim- 
ulations for the triangular lattice in the thermodynamic 
limit near the phase transition. The predictions from the 
effective model also show a sharp increase at the phase 
transition, but the agreement is limited to a smaller re- 
gion since the contribution of the spin current response 
was not included. In any case, the agreement near the 
phase transition line again gives strong support for the 
validity of the effective model in Eq. 



IV. ANALYSIS OF THE PHASE TRANSITION 

Finally, we would like to consider the implications of 
our model for the supersolid to superfluid transition, 
which has been reported to be second order before.- If 
our model is correct this transition should simply corre- 
spond to a melting of the solid due to the appearance 
of more and more virtual excitations, i.e. a breakdown 
of the perturbation expansion. While we have not found 
an energetic condition which would quantify the location 
of the supersolid to superfluid transition line, we believe 
that such a melting transition should be first order. One 
reason is also that a third order term in the Ginzburg 
Landau expansion is generically allowed by momentum 
conservation for an order parameter with triangular sym- 
metry, so that a continuous transition can only occur at 
special points. We have therefore tested the behavior of 
the phase transition with QMC simulations. In Fig. [5] it 
can be seen that for t = 0.16 the structure factor abruptly 
goes to zero close io p, ~ 2.012(2) while the densities 
show clear discontinuities. For larger p and smaller t 
closer to the particle-hole symmetric point these discon- 
tinuities become systematically smaller, but can still be 
observed even without histogram methods. The discon- 
tinuity in p corresponds to a region of phase separation 
if the density is held fixed as shown in the inset. During 
preparation of the paper we discussed our discovery of the 
first order behavior with one of the authors of Ref. Q and 
since then two quantitative studies have appeared which 
confirm the first order behaviori^ 



0.15 




0.05 



p-l/3 
Ps 

S(Q) 



13 0.14 0.15 0.16tA?0.17 



Supersolid 



Superfluid 



2.09 



2.1 



2.11 



FIG. 5: The density, structure factor S{Q) = 

(IEr'5're"^ ''P>/^^ at g =[4/37r,0] aird the superfluid den- 
sity with L = 36 and T = 0.01 at t = 0.16 and 1/ = 1 by 
varying /i. Inset: Observed region of phase separation be- 
tween supersolid and superfluid phases. 



can be calculated with the help of spin-wave theory for 
positive and negative hopping t. For positive t the calcu- 
lations agree well with QMC simulations. For t < we 
make a quantitative prediction for the supersolid region, 
which has so far not been estimated by numerical meth- 
ods. We find that the kinetic frustration for t < does 
not destroy the supersolid phase, but the higher order 
corrections imply a narrower region. In fact, the higher 
order terms imply a possible maximum of the extent of 
the supersolid phase with increasing \t\ in this case. It 
is indeed feasible that a configuration corresponding to 
a low density (spin down) on one sublattice and a small 
x-y antiferromagnetic order on the other two sublattices 
may be stable all the way to the SU(2) invariant point 
V = —2t = p/3. In this case the SU(2) invariant point 
would play the role of a trivial "degenerate" supersolid 
where a superposition of broken orders is possible in any 
direction. 
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V. CONCLUSIONS 

In conclusion we proposed an analytical model in the 
form of defects moving in an ordered background, which 
gives a quantitative description of the solid to supersolid 
transition in the hard-core boson model. In this way the 
transition line can be determined and correlations nearby 



Appendix: Spin wave analysis 

We consider a two dimensional system of hardcore 
bosons on an effective hexagonal lattice in Eq. (2). To en- 
force the hardcore constraint in a simple way, the Hamil- 
tonian is mapped onto the two dimensional XXZ model 



5 



with external magnetic field using aj O , Qi o S~ , 
and hi O 5*^^ + 1/2. With this mapping the hardcore 
boson model becomes 

H = -Y^UStS- + SrS+) + ^ y5f 5| - S J] 

i,d (ij) i 

(A.l) 

where B = — XV is an effective field and A = 3/2 
is half the coordination number on a hexagonal lattice. 
Here and in what follows we generally omit the ground 
state energy. 

The z-axis can be aligned with the mean field mag- 
netization direction by the rotation Sf = S'{^cos{9) + 
S'f sin(0), Sf = 5-^ and Sf = -Sf sin(6l) + cos(6l), 
where 5*^ is the spin vector in the rotated frame. The 
hexagonal lattice has two spins per unit cell, which 
are are expressed in terms of bosonic operators Ci (c|) 
and d, (4) respectively, i.e. 5f = ^{S'+ + S'^) = 
|(c| + Ci), Si^ = ^{Si'^ ~ Si ) = ^[{cl — Ci), S'i^ = 
^ — cjci and likewise for the other sublattice in terms 

of di (dl). Note that the spin- wave operators obey the 
usual (soft-core) bosonic commutation relations. Sub- 
stituting these expressions into Eq. (IA.1[) and assum- 
ing a dilute gas of spin- waves, we can ignore cubic 
and quartic terms in the bosonic operators. The lin- 
ear term in the operators is required to vanish yield- 
ing, cos6' = iJj- — XV) / {^XV + J2d^d) , where the sum 
runs over the nearest, next- nearest and next-next- nearest 
neighbor bonds. The longer range hopping terms of or- 
der play an important role in determining the mean 
field angle and the superfluid density below, but they can 
be neglected when determining the Bogoluibov rotation. 
Therefore, the Hamiltonian reduces to 



H = 



E 



{A[cld]+c,d,)+iA-t) (c 



cldj + Cid 



1)} 



+3t i clc + Y 4d^ ) + Hmf, 



(A.2) 



where A = sin^ Q (2t + V) j ^ and the mean- field Hamil- 
tonian is given by 



MF 



N 



cos^e 



cos 9 



-B 



(A.3) 



which is minimized by the same classical condition 
cos6' = B/ {XV -\-^j^td)- Here N is the total number 
of lattice sites. 

The Hamiltonian (|A.2p can be diagonalized by the 
Fourier transformation q — y^27iV e**'' '''Cfc and di = 
y/2/NJ2k e^'^'^'Ck- We introduce the structure factor 



7fe = 3 Ed 



'Ck- 
Jk.ra — 



"jklc^ which is a complex number. 



The phase (pk can be absorbed in a gauge transforma- 
tion dk — e~*'^'=(ife in the Fourier transformed Hamil- 
tonian and the amplitude of the structure factor reads 



\lk 



3 -t- 4 cos cos 



VSky 



2 cos (\/3fcy) 



•,1/2 



The 



Hamiltonian can be diagonalized using the Bogoluibov 
transformation^ 



Cfc = u^ckfe + Wfc alj. -I- uf /3fe 



4filk 



K^k + f^al^ - u^^/3k + v^pl^ 



where — cosh 61^, 



vl = sinh6'f , with tanh26'^ 



= sinhe^, 



— cosh 6*^, and 
lk\A/ (t-f |7fe|(A-t)) 
and tanh26'^ = -\-fk\A/ {t - \-ik\[A - t)). Finally, the 
Hamiltonian (|A.2p takes the diagonal form 



aloLk 



ll2)+Y.4{Plh + ll2) 



with cj^ 



6v/(i+|7/c|(v4-<))2-(A|7fe|)2 and = 



(A.4) 

&^{t~\^k\{A-t))^^{A\^k\)'. 

The boson density of the model in Eq. (2) is directly 
given by p = ((S*^) -I- 1/2) where the magnetization reads 

(5^) = (-5f sin 61 4- cos 61) 



cos 9 



which follows from a straight-forward evaluation of the 
expectation values in terms of the diagonal bosons in 
Eq. (TOl) . 

The superfiuid density is given by the second deriva- 
tive of the energy of the spin system with respect to a 
uniform twist v across the system ps = d"^ {H [v)) / d^ v . 
The Hamiltonian H{v) is derived from the XXZ model 
in Eq. ljA.ip by application of a site dependent ro- 
tation by an angle Vi around the z axis, Sf — 
S'+e*''-, S- SiC-^"^ and Sf Sf. Expand- 
ing the Hamiltonian around 
leads to H{v) = H + J2 

J! - 

and Td - (l/2)trf {S+Sr^, + S^ 
energy'" 

stiffness is given by ps = ^-§^Y.i,d ^l^^d)- Second or- 
der perturbation leads to a term integrating the current- 
current correlator with respect to Js, which can be ne- 
glected in our spin-wave approachi^S, To obtain the spin 
stiffness a uniform twist v along one direction is applied, 
which leads to 



- i^i+d = 
-iyj/2Td), where 

the spin current 
d ^ (i/z,)id + '^i+d^r) spin- kinetic 

To first-order in perturbation theory the spin 



i^d = 

^i+d ~ ^i+d^i ) 



-E 



^d ] ^d{SiSf^d + Sfj^j^Sf), 



OJ. 

'-"'3 / Qx Qx i ay Qy\ 



(A.5) 



The spin-spin correlation functions are provided by 

{SfS^ + Spy) = cos'' 9{S'fSf) + sin^9{S['Sf) + 
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{Sf' S'j'), where d = 1,2,3 denote nearest neighbor, next- and the functions fd{k) are given by 
nearest neighbor, and next-next-nearest neighbor, re- 
spectively. The spin-spin correlation functions in the ro- 
tated frame read fi{k) — 2|7fe|, 



(sr^;;,) = ((5'^)-i/4),vd, 



+2 cos{3k^/2) cos{3V3ky/2) 
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